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Ray stability is investigated in environments consisting of a range-independent background sound-speed pro- 
file on which a range-dependent perturbation is superimposed. Theoretical arguments suggest and numerical 
results confirm that ray stability is strongly influenced by the background sound speed profile. Ray instability 
is shown to increase with increasing magnitude of a(I) = (I / co)dco/dI , where 2n/co(I) is the range of a ray 
double loop and / is the ray action variable. This behavior is illustrated using internal-wave-induced scattering 
in deep ocean environments and rough surface scattering in upward refracting environments. 



o 
o 

(N 



Q 

U 
d 



(N 
> 
oo 
m 
o 

oo 
O 
CI 

o 



> 
X 



PACS numbers: 43.30.Cq, 43.30.Ft, 43.30.Pc 



I. INTRODUCTION 

Measurements made during the Slice89 propagation 
experiment 1 , made in the eastern North Pacific, showed a clear 
contrast between highly structured steep ray arrivals and rela- 
tively unstructured near-axial flat ray arrivals. Measurements 
made during the AET experiment^ in a similar environment 
provided further evidence of the same behavior. Motivated by 
these observations, several authors 4 - 5 - 6 have investigated ray 
sensitivity to environmental parameters. The work described 
in this paper continues the same line of investigation. 

Like the earlier work we focus on ray path stability in physi- 
cal space or phase space. The extension to travel time stability 
is not considered in this paper. Clearly, however, travel time 
stability must be addressed to fully understand the Slice89 and 
AET measurements. In spite of this limitation, and limita- 
tions of the ray approximation, our results are in agreement — 
qualitatively, at least — with the Slice89 and AET observa- 
tions. 

In Sec. [IT] we provide the theoretical background for the 
work that follows. Most of this material builds on the action- 
angle form of the ray equations. A trivial observation that fol- 
lows from the angle-action formalism is that ray sensitivity 
to the background sound speed structure is controlled by the 
function <X)(I), where 2n/co is the range of a ray double loop 
and / is the action variable. A heuristic argument suggests that 
dco/dl should be closely linked to ray stability. In Sec. ITTilwe 
present numerical simulations that are chosen to demonstrate 
the importance of the background sound speed structure on 
ray stability. Simulations are shown for both internal-wave- 
induced scattering in deep ocean environments and rough sur- 
face scattering in upward-refracting environments. The latter 
are included to demostrate the generality of the arguments pre- 
sented. These simulations strongly suggest that ray stability is 
controlled by the magnitude of the nondimensional quantity 
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In Sec. [IV] we explain the mechanism by which \a\ controls 



ray stability. In Sec. [VJwe summarize our results and briefly 
discuss: (i) the relationship between our work and earlier in- 
vestigations; (ii) timefront stability; (iii) the dynamical sys- 
tems viewpoint; and (iv) the extension to background sound 
speed structures with range-dependence. An appendix is re- 
served for some mathematical details. 



II. THEORETICAL BACKGROUND 

A. One-way ray equations 

We consider underwater sound propagation in a two- 
dimensional waveguide with Cartesian coordinates z (upward) 
and r (along-waveguide). One-way ray trajectories satisfy the 
canonical Hamilton's equations (cf. e.g. Ref. □ and refer- 
ences therein), 
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with Hamiltonian 



h(p,z;r) = -\ft 



-2 . 



Here, c(z;r) is the sound speed; p is the vertical slowness 
which is understood as the momentum, conjugate to the gener- 
alized coordinate z; and r is the independent (time-like) vari- 
able. The vertical slowness and the sound speed are related 
through pc = sin<p, where <p is the angle that a ray makes 
with the horizontal. 



B. Near-integrability under small waveguide perturbations 

Assume that the sound speed can be split as into a back- 
ground (range-independent) part, C(z), and a small range- 
dependent perturbation, 8c(z;r). Then, to lowest-order in 
8c/c, the Hamiltonian takes the additive form 

h = H(p,z) + 8h(p,z;r). 

Introduce now the Poincare action 



a ' Author to whom correspondence should be addressed. 

fberon@rsmas .miami . edu 



Electronic mail: 



— &dz P = - dzVc-t-H 2 , 
2n J n Jz_ 



(3) 



Typeset by REVTr^C 



2 



where z± denotes the vertical coordinate of the upper (+) and 
lower (— ) turning points of a ray double loop, and consider the 
canonical transformation into action-angle variables {p,z)>— > 
defined by 



c-^)-m. 



According to the above transformation, 

h(p,z;r) ^> H(I) + Sh(I,$;r) 
and the ray equations (0 take the form 

dl d ST dd d 



where 



, . dH 
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Set (0} constitutes a near-integrable Hamiltonian system for 
the reasons explained next. 

In the limit 8h — ► 0, the ray equations (0}, which have one 
degree of freedom, are autonomous and the corresponding 
Hamiltonian, H, is an integral of motion that constrains the 
dynamics. As a consequence, the equations are integrable 
through quadratures and the motion is periodic with (spatial) 
frequency (O. Namely I — Iq and ■& = ■&() + cor mod27T, where 
Io and #o are constants. Every solution curve is thus a line 
that winds around an invariant one-dimensional torus, whose 
representation in (p,z)-space is the closed curve given by the 
isoline H = H(Iq). Notice that C -1 cos <p = H; consequently, 
each torus can be uniquely identified by the ray axial angle de- 
fined by (p„(p,z) = arccosC„H, where C a is the background 
sound speed at the sound channel axis. 

With nonzero 8h, the Hamiltonian, H + 8h, is no longer 
an invariant (the equations are nonautonomous) and the sys- 
tem may be sensitive to initial conditions, leading to chaotic 
motion in phase space. The distinction between regular and 
chaotic ray trajectories is commonly quantified by the Lya- 
punov exponent, formally defined by 



1 d 
Voo = lim lim -In — , 

r->°°do->0 r do 



(6) 



where d(r), such that d(Q) — do, is a suitably chosen mea- 
sure of the separation between neighboring ray trajectories in 
phase space. For regular trajectories, d ~ r" as r — > °° with 
a > and, hence, Voo = 0. For chaotic trajectories, instead, 
d ~ exp Vc»r; in this case, v" 1 , the average e-folding range, is 
regarded as the predictability horizon. 



C. Variational equations 

Although any norm of (8p,Sz), a perturbation to a trajec- 
tory (p,z), could be used to define a distance in phase space, 
this is not a trivial task because z and p do not have the same 



dimensions, which complicates the evaluation of 0. The fol- 
lowing procedure eliminates this problem. 

The variational equations, which follow from the ray equa- 
tions (0, are 
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where 
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with po — p(Q) and zo = z(0), and I is the identity matrix. 
Here, J(r) and Q(r) are the Jacobian matrices of the Hamil- 
tonian vector field and associated flow, respectively; the last 
is usually referred to as the stability matrix. Notice that 
(8p,8z) T = Q(8po,8zo) T at the lowest-order in (8po,8zo). 
Let now v Q (r) be the largest of the two eigenvalues of Q, and 
consider the definition (cf. e.g. Ref. 1221) 
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If the limit in (|8|i exists, and is not nil, then (8p, 8z) ~ exp Vc*,r 
as r — > oa (nearby trajectories diverge exponentially in range) 
and, hence, (|8jl can be taken as a suitable definition of the Lya- 
punov exponent. The equivalence between © and (|8j can be 
understood by noting that the variational equations describe 
the evolution of an infinitesimal circle of initial conditions sur- 
rounding a specified ray initial condition. The circle gets de- 
formed into an ellipse whose area equals that of the initial cir- 
cle according to Liouville's theorem (cf. e.g. Refs. I8l9i> . The 
eigenvectors of Q define the orientation of the ellipse. The 
largest eigenvalue is a measure of the length of the semimajor 
axis and, hence, a suitable choice of d. 

A simple but very important observation follows from the 
action-angle formalism. Dependence of the ray and varia- 
tional equations on the background sound speed structure en- 
ters only through the function co(I). The action-angle form 
of the variational equations for the perturbed system strongly 
suggests that ray stability and dco/dl are closely linked. The 
mechanism through which dco/dl influences ray stability can 
be seen from the action-angle form of the ray variational 
equations, 
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If one assumes that the second derivatives of 8h are zero- 
mean random variables, then when do/d/ = these terms 
should lead to slow (power-law) growth of 8~& and 81. But 
if |do/d/| is large, this term will cause \8~&\ to rapidly grow 
for any nonzero 1 8l\ . The perturbation terms will then lead to 
a mixing of and |5/|. The term dco/dl will lead, in turn, 
to further growth of |5t?| . As this process repeats itself, both 
|5/| and |5t?| are expected to grow rapidly. Thus ray instabil- 
ity is expected to be significantly enhanced when |da)/d/| is 
large. 
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D. KAM theory 

A (background) Hamiltonian is said to be nondegenerate 
if dco/dl 7^ 0, i.e. if the frequency varies from torus to torus, 
which is a condition for nonlinearity of the system. An im- 
portant result for near-integrable nondegenerate Hamiltonian 
systems is the celebrated Kolmogorov-Arnold-Moser (KAM) 
theorem on the stability of periodic solutions; cf. e.g. Refs. 
I8l9l This theorem states that if 8h is small enough, for most 
initial conditions the motion remains periodic (i.e. confined 
to tori) and the complement of the periodic motion (i.e. the 
chaotic motion) has a measure that tends to zero as 8h — ► 0. 
The KAM theorem thus guarantees that periodic motion (i.e. 
the KAM tori) separate the destroyed tori, leading to the no- 
tion of "islands" of (eternal) stability immersed in a chaotic 
"sea." 

The mechanism that produces the destruct ion of tori is 
trajectory-medium resonance; cf. e.g. Refs. I8l9i For ex- 
ample, assume the perturbation Hamiltonian 8h to be peri- 
odic in range with frequency £1. Then it can be represented 
in Fourier series 8h = Re£ m „A m „(/)expi(ra# —n£lr). For 
KAM or nonresonant tori, mco + nQ. ^ for all integers n,m 
and the motion is periodic. In contrast, those tori that sat- 
isfy mco +n£l — for some integers n,m are said to be in 
resonance and chaotic motion develops. If several tori are 
captured into resonance, then the character of the chaotic 
motion will depend on whether these resonances overlap or 
not. For instance, consider two resonances centered at I\ 
and h- The widths of these resonances can be estimated as 
Alj = 4a/A// |da),-/d/|, where A; is the amplitude of the reso- 
nant term in the above expansion. Define then A/ = A/i + A/2. 
A benign form of chaos is present when these resonances are 
isolated, i.e. AI < \I\ — h\', strong chaos, in contrast, emerges 
when these resonances overlap, i.e. AI > \I\ — The last 
criterion, due to Chirikov 11 , gives thus a quantitative estimate 
of the size of the region of phase space occupied by chaotic 
trajectories. 

The focus in KAM theory on the role of individual ray- 
medium resonances might seem to be at odds conceptually 
with the heuristic argument given at the end of the preceding 
subsection. There it was argued that the second derivatives of 
the perturbation to the environment could be treated as ran- 
dom variables. Ref. [3 partially bridges this conceptual gap 
by showing that the KAM theorem can be applied to prob- 
lems for which the perturbation consists of a superposition of 
an arbitrarily large finite number of frequencies. Chirikov's 
criterion is still applicable, but its evaluation seems feasible 
only if the number of frequencies that comprises the perturba- 
tion is very small. 



III. NUMERICAL SIMULATIONS 

In this section numerical simulations are presented which 
were chosen to demonstrate the importance of the background 
sound speed structure on ray stability. We consider first deep- 
water conditions and upward-refracting conditions afterward. 




FIG. 1: Background sound speed profiles used in the numerical work 
presented in this paper. 



A. Deep-water conditions 

Four different background profiles are studied here (cf. Fig. 
G}. These are designated S89, C89, C18, and cosh. The S89 
profile corresponds to a range average of the Slice89 sound 
speed observations. The C89 profile is a canonical profile 13 
with parameters chosen to approximately match S89's sound 
channel axis depth, axial sound speed, and surface sound 
speed. The C18 profile can be regarded as an idealized model 
of the sound speed structure in the North Atlantic, whose up- 
per ocean structure is associated with the 18 °C water mass 7 . 
Finally, the cosh profile, with cosh-dependence on depth rela- 
tive to the axis, was chosen because it has the special property 
da)/d/ = for all /. 

The same internal-wave-induced sound speed perturbation 
field is superimposed on all four background profiles. This 
field is assumed to satisfy the relationship 

8c/c = iiN 2 C, (9) 

where £(z;r) is the internal-wave-induced vertical displace- 
ment of a water parcel and N(z',r) is the Brunt-Vaisala fre- 
quency. Relationship (|9} with jU = 1 .25 s 2 m~ 1 was found to 
give a good fit to AET hydrographic measurements^ Our 
simulated internal-wave-induced sound speed perturbations 
are similar to those used in Refs. these are based on 

Eq. l|9} and make use of the N profile estimated from mea- 
surements made during the AET experiment. The statistics 
of £ are assumed to be described by the empirical Garrett- 
Munk spectrum '4. The vertical displacement £ is computed 
using Eq. (19) of Ref. [H with the variable x replaced by r 
and y = = t. Physically this corresponds to a frozen verti- 
cal slice of the internal wave field that includes the influence 
of transversely propagating internal wave modes. A Fourier 
method is used to numerically generate the sound speed per- 
turbation fields. A mode number cutoff of 30 and a horizontal 
wavenumber cutoff of 2n km were used in our simulations. 
The internal wave strength parameter was taken to be the nom- 
inal Garrett-Munk value. 

Fig. |2 shows the stability exponent, v (a finite-range esti- 
mate of the Lyapunov exponent; cf. appendix A), as a function 
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FIG. 2: Stability exponents (finite range estimates of Lyapunov ex- 
ponents) as a function of the initial depth and launch angle, for ray 
simulations in the (a) C89, (b) CI 8, (c) S89, and (d) cosh waveguides 
(cf. Fig. 1) with internal- wave-induced perturbations superimposed. 



of initial ray position in phase space for each of the four envi- 
ronments considered. This figure provides a general picture of 
the ray motion stability character in each of the waveguides. 
In the C89 waveguide, ray trajectories with small (resp., large) 
unperturbed ray axial angles have small (resp., large) asso- 
ciated stability exponents. This trend is reversed in the S89 
waveguide. The disparity in the stability properties of the ray 
motion in these waveguides contrasts with the close similar- 
ity of the corresponding background sound speed profiles. In 
the C18 waveguide ray trajectories have in general relatively 
small associated stability exponents, except in a narrow band 
of initial actions (or unperturbed axial angles) where the ex- 
ponents are large and the ray motion more unstable. In oppo- 
sition to the other waveguides, in the cosh profile the stability 
exponents are very small everywhere in phase space and the 
ray motion is remarkably stable. 

In each panel of Fig. [2] the horizontal line shown corre- 
sponds to a fan of rays, launched on the sound channel axis, 
with positive angles. For these rays, in each environment ray 
depth at a range of 1000 km and stability exponent are plot- 
ted as a function of launch angle, <p , in Fig. [5] In that fig- 
ure bands of regular trajectories appear as a smooth slowly- 
varying curve z(<Pq) with small associated values of v; bands 
of chaotic rays appear as a highly structured z(<p ) curve with 
large associated values of v. Also shown in Fig. [3]is a plot of 
the stability parameter a defined in Eq. Q vs. <p in each en- 
vironment. (The action / is a monotonically increasing func- 
tion of the axial ray angle which, for rays starting at the axis, 
coincides with <p ; thus replacing / by <p represents a simple 
stretching of the abscissa.) Notice in the z vs. <p plots the 
seemingly unstructured (resp., ordered) distribution of points 
associated with those angular bands where v is large (resp., 
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FIG. 3: Each panel shows ray final depth (upper-plot), stability expo- 
nent (middle-plot), and stability parameter (lower-plot) at 1000 km 
range as a function of ray launch angle for a source on the sound 
channel axis. The sound speed structures are identical to those used 
in the panels of Fig. 2 with the same labels. 



small). Notice also that v (or, equivalently, the irregularity of 
Z as a function of <p ) appears to increase with increasing \ a\. 
In particular, note that a = (since dco /dl ~ 0) for all <p in 
the cosh profile. In this case ray final depth varies very slowly 
with <p , and the stability exponents are very small. 

The numerical results presented in Figs. [5] and [3] suggest 
that ray stability is strongly influenced by the background 
sound speed structure and that ray instability increases with 
increasing \a \ . These results are consistent with the heuristic 
argument given a the end of Sec. Qjjbased on the action-angle 
form of the ray variational equations; there it was argued that 
ray instability should increase with increasing |do/d/| . 



B. Upward-refracting conditions 

The validity of the argument given in Sec. [H]is not limited 
to deep ocean conditions, strongly suggesting that the result 
should be more generally valid. We now present numerical 
results that support this expectation. 

Figure[4]shows two upward-refracting sound speed profiles. 
Figure |3 shows a in the same environments and the differ- 
ence in range, Ar, between perturbed (rough surface) and un- 
perturbed (flat surface) rays as a function of launch angle at 
the surface after making 21 loops, which corresponds to 20 
surface reflections. [In this type of environment the defini- 
tion of / in Eq. l|3} is unchanged except that the upper in- 
tegration limit is z = for all rays.] Rough surface scattering 
was treated using a frozen simulated surface gravity wavefield 
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FIG. 4: Background sound speed profiles used to construct the curves 
in Fig. 5. 

with a kr 1 ! 1 surface elevation wavenumber spectrum with 
0.02 rad trr 1 < k < 0.16 rad m _1 , Ak = 10~ 3 rad m -1 , and 
rms slope of 4 x 10~ 3 . To treat specular ray reflections from 
this surface, the surface boundary condition was linearized; 
the surface elevation was neglected, but the nonzero slope was 
not approximated. 

Figure[5]shows clearly that ray stability, as measured by Ar, 
is controlled almost entirely by the background sound speed 
structure via a, rather than details of the rough surface. 



IV. SHEAR-INDUCED RAY INSTABILITY 

In this section we present additional numerical simulations 
that give insight into the mechanism by which a influences 
ray instability. For simplicity we assume deep water condi- 
tions here. The arguments presented here make use of the 
well-known (cf. e.g. Refs. 1819 ) analogy between ray motion, 

I Vol Ivo 
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FIG. 5: Upper panels: absolute value of the difference between per- 
turbed (rough surface) and unperturbed (flat surface) ray range after 
20 surface reflections as a function of initial ray angle at the surface 
for sound speed profiles PA (a) and PB (b) in Fig. 4. Lower pan- 
els: stability parameter as a function of ray angle in each of the two 
environments. 



defined by Eqs. (0 and particle motion in an incompressible 
two-dimensional fluid. 

The left panels of Fig. |6| show the range evolution of a 
segment of a Lagrangian manifold (a smooth curve in phase 
space) in the unperturbed C89 (top) and cosh (bottom) waveg- 
uides. The segment is depicted in black at r = and in red at 
r = 1000 km in both waveguides. As a consequence of Li- 
ouville's theorem the segment cannot break or intersect itself 
but it can increase in complexity as range increases. In the un- 
perturbed case, since the motion is integrable (i.e. each point 
of the segment preserves its initial /), the length of the seg- 
ment can grow in range, at most, following a power law. The 
cosh profile has a special property. In that profile the mani- 
fold segment just rotates counterclockwise at a constant rate 
co. The reason for the difference in behavior is that in the 
C89 waveguide CO varies with /, whereas in the cosh profile co 
does not vary with /. The monotonic decay of CO as a function 
of / in the C89 waveguide induces a "shear" in phase space 
which causes the outermost points of the segment to rotate 
more slowly than the innermost ones and, hence, causes the 
segment to spiral. The ray motion in phase space associated 
with the unperturbed C89 waveguide can thus be regraded as 
analogous to that of ideal fluid particles passively advected by 
a stationary planar circular flow with radial shear. In the cosh 
profile there is no shear. In polar coordinates radial shear can 
be defined as 

where p is the radial coordinate and uq is the -component 
of the velocity field. (More correctly, this quantity is, apart 
from a factor of 2, the p9 -component of the strain-rate tensor 
for planar circular flow; cf. e.g. Ref. [3) The connection 
with ray motion in phase space can be accomplished by iden- 
tifying / with p and col with uq. The replacements p i — >■ 7 
and uq col in Eq. dlOt thus give the analogous expression 
Idco/dl for the shear in phase space. Notice that this expres- 
sion is (apart from the co -1 -factor) the stability parameter a. 
We have chosen to include the o -1 -factor in the definition of 
a because of precedent 10 and because it is convenient to make 
a dimensionless. 

The right panels of Fig. [6] show the evolution of the La- 
grangian manifold segment in the same waveguides as those 
used to produce the left panels but with a superimposed per- 
turbation induced by internal wave fluctuations. Notice the 
highly complicated structure of the Lagrangian manifold in 
the perturbed C89 waveguide as compared to that in the un- 
perturbed one. (Note that the fan of rays used to produce the 
figure is far too sparse to resolve what should be an unbro- 
ken smooth curve which does not intersect itself.) In contrast, 
observe that in the cosh environment the sound speed pertur- 
bation has only a very minor effect on the evolution of the 
Lagrangian manifold. 

Perturbations to steeper rays caused by internal-wave- 
induced sound speed perturbations in deep ocean, including 
those in our simulations, are significant only near the ray's 
upper turning depth. This observation motivates a simple 
model that gives insight into the mechanism by which a is 
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FIG. 6: Evolution of a portion of Lagrangian manifold in the C89 
(upper panels) and cosh (lower panels) waveguides with (right pan- 
els) and without (left panels) internal- wave-induced perturbations. 
Dashed and solid/dotted curves show the manifold at r = and 
r = 1000 km, respectively. 



linked to ray stability. In the model, each portion of a La- 
grangian manifold acquires a sinusoidal "wrinkle" at each up- 
per turning point, but is otherwise unaffected by the sound 
speed perturbations. The evolution of small segments of a La- 
grangian manifold using such a model, in both the C89 and 
cosh waveguides, is shown in Fig. Notice the rapid growth 
in complexity of the segment as the range increases in the C89 
waveguide. After acquiring a wrinkle, the segment stretches 
and folds as a result of the radial shear in phase space. As 
additional wrinkles are acquired this process is repeated suc- 
cessively in range, making the shape of the Lagrangian man- 
ifold segment even more complex. In opposition to this sit- 
uation, observe the simplicity of the segment's shape in the 
cosh waveguide as range increases. In this case, after acquir- 
ing a wrinkle, the distorted segment rotates counterclockwise, 
without the additional influence of shear-induced stretching, 
at a constant frequency co. 

Each time a perturbation is introduced, the action / changes 
by the amount 81, say, which we assume to be of the same 
order as the perturbation. As a consequence, to lowest-order 
R(= 2%/ 03) experiences the change 
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FIG. 7: Evolution of a segment of Lagrangian manifold in the C89 (a) 
and cosh (b) waveguides under the influence of an idealized range- 
dependent perturbation. The perturbations are is assumed to be in 
the form of a series of "kicks" that produce a sinusoidal deformation 
to the Lagrangian manifold at each upper turning depth. The figure 
shows a portion of Lagrangian manifold, originally located on a torus 
with frequency co(= 2k /R), immediately after experiencing a kick at 
range r=(k— l)R, k= 1,- •• ,4, and before receiving the next one at 
r = IcR, as well as at several intermediate stages. 



(l-a8l/I)R. 



The perturbation to the range of a ray double loop R depends 
on both the perturbation and the properties of the background 
sound speed structure. Under the change 1 1— > I + 5/, a suf- 
ficient condition for R to remain invariant at lowest-order is 
a = 0. This provides an explanation for the remarkable stabil- 
ity of the cosh waveguide. That is, when a = the ray motion 
remains periodic at lowest-order — no matter the complexity 
of the perturbation mounted on the waveguide. To lowest- 
order, a nonvanishing shear (a/0) appears as a necessary 
condition to sustain the successive stretching and folding of 



the Lagrangian manifold after it gets distorted by the pertur- 
bation. (Of course, chaotic motion is still possible when a = 
provided that the perturbation strength is sufficiently large.) It 
is thus expected that where \a\ is small (resp., large) there 
will be less (resp., more) sensitivity to initial conditions and, 
hence, the motion be more regular (resp., chaotic). Support 
for this conjecture is given in the numerical simulations pre- 
sented in this paper. 
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V. SUMMARY AND DISCUSSION 

In this paper we have considered ray motion in environ- 
ments consisting of a range-independent background sound 
speed profile on which a weak range-dependent perturbation 
is superimposed. The results presented show that ray sta- 
bility is strongly influenced by the background sound speed 
structure; ray instability was shown to increase with increas- 
ing magnitude of a (I) = (//£»)dO)/d/, where 2k/(o{I) is the 
range of a ray double loop and / is the ray action variable. 
This conclusion is based largely on numerical simulations, but 
the simulations were shown to support a simple heuristic ar- 
gument based on the action-angle form of the ray variational 
equations. The mechanism by which a controls ray instability 
was shown to be shear-induced enhancement of perturbations 
caused by the sound speed perturbation term. The importance 
of a was illustrated with numerical simulations of ray motion 
in deep ocean environments including internal-wave-induced 
scattering, and in upward-refracting environments including 
rough surface scattering. So far as we are aware, this conclu- 
sion is consistent with all of the numerical results presented 
earlier 4 ' 5,6 . Ref. 6 relied heavily on results that follow from a 
dynamical systems viewpoint. 

The connection between our work and results relating to 
dynamical systems deserves further comment. Recall that 
the condition do/d/ 7^ (the nondegeneracy condition) must 
be satisfied for the KAM theorem to apply, and that this re- 
sult guarantees that some rays are nonchaotic provided the 
strength of the range-dependent perturbation is sufficiently 
weak. This theorem might seem to conflict with our assertion 
that ray instability increases with increasing \a\. This appar- 
ent conflict can be resolved by interpreting our statement as 
a statement of what happens for most rays. That is, for most 
rays stability exponents (finite range estimates of Lyapunov 
exponents) increase as \a\ increases. This does not rule out 
the possibility that for a fixed dco/dl ^ some rays will be 
nonchaotic. 

A seemingly more troublesome conflict between our sim- 
ulations and KAM theory follows from the result, noted ear- 
lier, that each isolated resonance has a width proportional to 
|d<a/d/| -1 / 2 . Doesn't this imply that rays should become in- 
creasingly chaotic as |do/d/| decreases? The answer, we 
believe, is no. To understand why, consider rays in a band 
^0 < I < h over which da)/d/ is approximately constant, 
dco/dl = (<Bi — <Bo)/(^i — h)- Within this band resonances 
are excited at selected values of co. The number of resonances 
excited is approximately proportional to — ©o|, which, 
in turn, is proportional to |dto/d/|. The "degree of chaos" 
should be proportional to the product of the number of reso- 
nances excited and the width of individual resonances. This 
product scales like |d©/d/| 1//2 ; this suggests that rays should 
become increasingly chaotic as |d<»/d/| increases, consistent 
with our simulations. Neither the KAM theorem nor the reso- 
nance width estimate estimate applies in the limit dco/dl = 0. 
Behav ior in that limit is likely problem-dependent (e.g. Refs. 
Il9l20l) . Our simulations — which are probably representative 
of problems characterized by an inhomogeneous background 
and a weak perturbation with a broad spectrum — suggest that 



ray motion in this limit is very stable. 

In this paper we have addressed the issue of ray stabil- 
ity in physical space or phase space; we have not addressed 
the related problem of travel time stability. The latter prob- 
lem is more difficult inasmuch as that problem involves, in 
addition to the (one-way) ray equations Q, a third equa- 
tion dT/dr = L and imposition of an eigenray constraint. 
Here, L — pdz/dr — h(p,z',r) using standard variables, or 
L = Idd-/dr — h(I,Td-;r) using action-angle variables. We have 
seen some numerical evidence that ray instability in phase 
space is linked to large time spreads, but this connection is 
currently not fully understood. 

An advantage of our use of the action-angle formalism 
is that essentially the same results apply if the assumption 
of a background range-independent sound speed profile, i.e. 
C = C(z), is relaxed to allow for a slowly-varying (in range) 
background sound speed structure, i.e. C = C(z',£r). (Here, e 
scales like the ratio of a typical correlation length of the range- 
dependence to a typical ray double loop range.) In the latter 
case the action is not an exact ray invariant but it is an adia- 
batic invariant, i.e. dl /dr = 0(e 2 ). Thus, correct to 0(e) the 
latter problem can be treated as being identical to the former 
one. Consequently, the problem that we have treated here also 
applies to slowly-varying background environments. 
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APPENDIX: STABILITY EXPONENTS 

Liouville's theorem (cf. e.g. Refs. I8l9l) guarantees that 
areas in the present two-dimensional phase space are pre- 
served, leading to detQ = 1 as a corollary. Consequently, 
2v± = trace Q ± \Arace 2 Q- 1. (It can be shown that v± — 
trace* Q as r — > °o; accordingly, Eq. (jSJl can be replaced 
by the equivalent expression v = Hindoo In |traceQ| as 
in Ref. 20.) The condition detQ = 1, however, is difficult to 
fulfill in highly chaotic flows due to the limitation of machine 
numerical precision. More precisely, an initial area in phase 
space tends to align along the one-dimensional perturbation 
subspace spanned by the eigenvector of Q associated with 
the largest Lyapunov exponent, making the elements of Q ill- 
conditioned. This implies that neither det Q nor trace Q — and, 
hence, neither v Q nor Voo — can be computed reliably at long 
range. The common approach to overcome this problem in- 
volves successive normalizations of the elements of Q, while 
integrating simultaneously Q and 0, after a fixed number of 
range steps 22 . An alternative strategy, which does not involve 
normalizations, has been taken here. 
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Let A(r) be the bounded-element matrix such that 

Q = Ae^, (11) 

which implies 

y Q = y A e Ar 

Here, A is a guessed value of Vo». The reason for introducing 
the decomposition (lilt is that if A is close to Voo, the matrix 
A will remain well conditioned at ranges long past those at 



which Q becomes poorly conditioned. This leads, in turn, to 
significantly improved numerical stability. Eq. Q leads to the 
modified variational equations 

dA 

— = (J-AI)A, A(0) = l. (12) 

Notice that detQ = detAexpAr = 1 and, hence, detA ~ as 
r — > oo. Eq. d!2t has been integrated in this paper after choos- 
ing A, in order to compute a finite range estimate of Vc*>, which 
we call a stability exponent and denote by v. 
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